package math2;

/**
 * Gauss points on triangle refer from Gid
 * 
 * @author hbui
 * 
 */
public class GaussTriangleGid {

	private static final double tri3a = 0.09157621;
	private static final double tri3b = 0.81684757;
	private static final double tri3c = 0.44594849;
	private static final double tri3d = 0.10810301;
	private static final double tri3w1 = 0.054975871827661;
	private static final double tri3w2 = 0.1116907948390055;

	private static final double[][] gp5 = new double[][] {
			{ -0.333333333333333, -0.333333333333333, 0.450000000000000 },
			{ -0.059715871789770, -0.059715871789770, 0.264788305577012 },
			{ -0.059715871789770, -0.880568256420460, 0.264788305577012 },
			{ -0.880568256420460, -0.059715871789770, 0.264788305577012 },
			{ -0.797426985353088, -0.797426985353088, 0.251878361089654 },
			{ -0.797426985353088, 0.594853970706174, 0.251878361089654 },
			{ 0.594853970706174, -0.797426985353088, 0.251878361089654 } };
	private static final double[][] gp6 = new double[][] {
			{ -0.501426509658180, -0.501426509658180, 0.233572551452758 },
			{ -0.501426509658180, 0.002853019316358, 0.233572551452758 },
			{ 0.002853019316358, -0.501426509658180, 0.233572551452758 },
			{ -0.873821971016996, -0.873821971016996, 0.101689812740414 },
			{ -0.873821971016996, 0.747643942033992, 0.101689812740414 },
			{ 0.747643942033992, -0.873821971016996, 0.101689812740414 },
			{ -0.379295097932432, 0.273004998242798, 0.165702151236748 },
			{ 0.273004998242798, -0.893709900310366, 0.165702151236748 },
			{ -0.893709900310366, -0.379295097932432, 0.165702151236748 },
			{ -0.379295097932432, -0.893709900310366, 0.165702151236748 },
			{ 0.273004998242798, -0.379295097932432, 0.165702151236748 },
			{ -0.893709900310366, 0.273004998242798, 0.165702151236748 } };

	public static final double[][][] GAUSS_PTS = { //
			{ { 1.0 / 3, 1.0 / 3, 1.00 / 2 } }, //
			{ { 0.5, 0.0, 1.00 / 6 }, { 0.5, 0.5, 1.00 / 6 }, { 0.0, 0.5, 1.00 / 6 } }, //
			{ { tri3a, tri3a, tri3w1 }, { tri3b, tri3a, tri3w1 }, { tri3a, tri3b, tri3w1 },
					{ tri3c, tri3d, tri3w2 }, { tri3c, tri3c, tri3w2 }, { tri3d, tri3c, tri3w2 } }, //
			// below here is not supported in Gid
			{ { (gp5[0][0] + 1) / 2, (gp5[0][1] + 1) / 2, gp5[0][2] / 4 },
					{ (gp5[1][0] + 1) / 2, (gp5[1][1] + 1) / 2, gp5[1][2] / 4 },
					{ (gp5[2][0] + 1) / 2, (gp5[2][1] + 1) / 2, gp5[2][2] / 4 },
					{ (gp5[3][0] + 1) / 2, (gp5[3][1] + 1) / 2, gp5[3][2] / 4 },
					{ (gp5[4][0] + 1) / 2, (gp5[4][1] + 1) / 2, gp5[4][2] / 4 },
					{ (gp5[5][0] + 1) / 2, (gp5[5][1] + 1) / 2, gp5[5][2] / 4 },
					{ (gp5[6][0] + 1) / 2, (gp5[6][1] + 1) / 2, gp5[6][2] / 4 } },//
			{ { (gp6[0][0] + 1) / 2, (gp6[0][1] + 1) / 2, gp6[0][2] / 4 },
					{ (gp6[1][0] + 1) / 2, (gp6[1][1] + 1) / 2, gp6[1][2] / 4 },
					{ (gp6[2][0] + 1) / 2, (gp6[2][1] + 1) / 2, gp6[2][2] / 4 },
					{ (gp6[3][0] + 1) / 2, (gp6[3][1] + 1) / 2, gp6[3][2] / 4 },
					{ (gp6[4][0] + 1) / 2, (gp6[4][1] + 1) / 2, gp6[4][2] / 4 },
					{ (gp6[5][0] + 1) / 2, (gp6[5][1] + 1) / 2, gp6[5][2] / 4 },
					{ (gp6[6][0] + 1) / 2, (gp6[6][1] + 1) / 2, gp6[6][2] / 4 },
					{ (gp6[7][0] + 1) / 2, (gp6[7][1] + 1) / 2, gp6[7][2] / 4 },
					{ (gp6[8][0] + 1) / 2, (gp6[8][1] + 1) / 2, gp6[8][2] / 4 },
					{ (gp6[9][0] + 1) / 2, (gp6[9][1] + 1) / 2, gp6[9][2] / 4 },
					{ (gp6[10][0] + 1) / 2, (gp6[10][1] + 1) / 2, gp6[10][2] / 4 },
					{ (gp6[11][0] + 1) / 2, (gp6[11][1] + 1) / 2, gp6[11][2] / 4 } }
	//
	};

	/**
	 * get number of Gauss points for a specific order
	 * 
	 * @param p
	 *            order
	 * @return number of Gauss points associated with order i
	 */
	public static int getNumGaussPoints(int p) {
		if (p > GAUSS_PTS.length) {
			throw new Error(
					"exceed maximum order of Gauss integration rule on triangle. The maximum order is "
							+ p + " only");
		}
		return GAUSS_PTS[p - 1].length;
	}

	/**
	 * get xi1 coordinate of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return xi1
	 */
	public static double getGaussXi1(int p, int i) {
		return GAUSS_PTS[p - 1][i][0];
	}

	/**
	 * get xi2 coordinate of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return xi2
	 */
	public static double getGaussXi2(int p, int i) {
		return GAUSS_PTS[p - 1][i][1];
	}

	/**
	 * get weight of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return weight
	 */
	public static double getGaussWeight(int p, int i) {
		return GAUSS_PTS[p - 1][i][2];
	}

}
